Clustering of Aerosols in Atmospheric Turbulent Flow 
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A mechanism of formation of small-scale inhomogeneities in spatial distributions of aerosols and 
droplets associated with clustering instability in the atmospheric turbulent flow is discussed. The 
particle clustering is a consequence of a spontaneous breakdown of their homogeneous space distri- 
bution due to the clustering instability, and is caused by a combined effect of the particle inertia 
and a finite correlation time of the turbulent velocity field. In this paper a theoretical approach 
proposed in Phys. Rev. E 66, 036302 (2002) is further developed and applied to investigate the 
mechanisms of formation of small-scale aerosol inhomogeneities in the atmospheric turbulent flow. 
The theory of the particle clustering instability is extended to the case when the particle Stokes 
time is larger than the Kolmogorov time scale, but is much smaller than the correlation time at the 
integral scale of turbulence. We determined the criterion of the clustering instability for the Stokes 
number larger than 1. We discussed applications of the analyzed effects to the dynamics of aerosols 
and droplets in the atmospheric turbulent flow. 
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I. INTRODUCTION 

It is known that turbulence enhances mixing (see, e.g., 
[1-9]). However, numerical simulations, laboratory ex- 
periments and observations in the atmospheric turbu- 
lence revealed formation of long-living inhomogeneities 
in spatial distribution of aerosols and droplets in turbu- 
lent fluid flows (see, e.g., [10-21]). The origin of these 
inhomogeneities is not always clear but their influence 
on the mixing can be hardly overestimated. 

It is hypothesized that the atmospheric turbulence en- 
hances the rate of droplet colHsions (see, e.g., [18, 22- 
25]). In particular, the turbulence causes formation of 
small-scale droplet inhomogeneities, and it also increases 
the relative droplet velocity. In addition, the turbulence 
affects the hydrodynamic droplet interaction. The latter 
increases the rate of droplet collisions. These effects are 
of a great importance for understanding of rain forma- 
tion in atmospheric clouds. In particular, these effects 
can cause the droplet spectrum broadening and acceler- 
ation of raindrop formation [18, 23]. Note that clouds 
are known as zones of enhanced turbulence. The pref- 
erential concentration of inertial particles (particle clus- 
tering) was recently studied in numerical simulations in 
[26-28]. The formation of network- like regions of high 
particle number density was found in [28] in high resolu- 
tion direct numerical simulations of inertial particles in 
a two-dimensional turbulence. 

The goal of this study is to analyze the particle-fluid 
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interaction leading to the formation of strong inhomo- 
geneities of aerosol distribution due to a particle clus- 
tering instability. The particle clustering instability is a 
consequence of a spontaneous breakdown of their homo- 
geneous space distribution. As a result, at the nonlinear 
stage of the clustering instability, the local density of 
aerosols may rise by orders of magnitude and strongly 
increase the probability of particle-particle collisions. 

It was suggested in [29-31] that the main cause of the 
particle clustering instability is their inertia: the particles 
inside the turbulent eddies are carried out to the bound- 
ary regions between the eddies by the inertial forces. 
This mechanism of the preferential concentration acts in 
all scales of turbulence, increasing toward small scales. 
Later, this was contested in [32, 33] using the so-called 
"Kraichnan model" [34] of turbulent advection by the 
delta-correlated in time random velocity field, whereby 
the clustering instability did not occur. 

However, it was shown in [35] that accounting for a 
finite correlation time of the fluid velocity field results in 
the clustering instability of inertial particles. Note that 
the particle inertia results in the compressibility of par- 
ticle velocity field. The effects of compressibility of the 
velocity field on formation of small-scale inhomogeneities 
in spatial distribution of particles were first discussed in 
[36, 37]. In this study a theoretical approach proposed 
in [35] is further developed and applied to investigate the 
mechanisms of formation of small-scale aerosol inhomo- 
geneities in the atmospheric turbulent flow. In particular, 
we extended the theory of particle clustering instability 
to the case when the particle Stokes time is larger than 
the Kolmogorov time scale, but is much smaller than the 
correlation time at the integral scale of turbulence. 

Remarkably, the particle inertia also results in forma- 
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tion of the large-scale inliomogcneitics in the vicinity of 
the temperature inversion layers due to excitation of the 
large-scale instability (see [29, 32, 38]). This effect is 
caused by additional non-diffusive turbulent flux of par- 
ticles in the vicinity of the temperature inversion (phe- 
nomenon of turbulent thermal diffusion). The charac- 
teristic time of excitation of the large-scale instability of 
concentration distribution of aerosols varies in the range 
from 0.3 to 3 hours depending on the particle size and 
parameters of the atmospheric turbulent boundary layer 
and the temperature inversion layer. The phenomenon 
of turbulent thermal diffusion was recently detected ex- 
perimentally using two very different turbulent flows cre- 
ated by oscillating grids turbulence generator [39-41] and 
multi-fan turbulence generator [42] for stably and unsta- 
bly stratified fiuid flows. 

The paper is organized as follows. In Sec. II we present 
governing equations and a qualitative analysis of the clus- 
tering instability that causes formation of particle clus- 
ters in a turbulent flow. In Sec. Ill we estimate the scal- 
ings of the particle velocity in the turbulent fluid for the 
case when the particle Stokes time is much larger than 
the Kolmogorov time scale, but is much smaller than the 
correlation time at the integral scale of turbulence. In 
Sec. IV we perform a quantitative analysis for the clus- 
tering instability of the second moment of particle num- 
ber density for St > 1, where St is the Stokes number. 
This allows us to generalize the criterion of the cluster- 
ing instability obtained in [35]. Finally, in Sec. V we 
overview the nonlinear effects which lead to saturation of 
the clustering instability and determine the particle num- 
ber density in the cluster. In Sec. V we perform numeri- 
cal estimates for the dynamics of aerosols and droplets in 
atmospheric turbulent flow. The conclusions are drawn 
in Sec. VI. The detail analysis of the scalings of the par- 
ticle velocity in the turbulent fluid is given in Appendix 
A. The detail analysis of the clustering instability of the 
inertial particles is given in Appendix B. 



II. GOVERNING EQUATIONS AND 
QUALITATIVE ANALYSIS OF PARTICLE 
CLUSTERING 

To analyze dynamics of particles we use the standard 

continuous media approximation, introducing the num- 
ber density field n(t, r) of spherical particles with radius 
a. The particles are advected by an incompressible tur- 
bulent velocity field u{t, r). The particle material density 
Pp is much larger than the density p of the ambient fluid. 
For inertial particles their velocity v(f, r) =/: u(i, r) due 
to the particle inertia and divv(t,r) ^ (see [43, 44]). 
Therefore, the compressibility of the particle velocity 
field v(f,r) must be taken into account. The growth 
rate of the clustering instability, 7, is proportional to 
(|divv(i,r)|2) (see [29, 30, 37]), where (•) denotes en- 
semble average. 

Let 0(t,r) be the deviation of the instantaneous par- 



ticle number density n(t,r) from its uniform mean value 
= (n): e(t, r) = n(t, r) - A^. The pair cor- 
relation function of Q(t,r) is defined as $(R, r,t) = 
(8(t, r -|- R)0(t, r)). For the sake of simplicity we will 
consider only a spatially homogeneous, isotropic case 
when <&(R, r, t) depends only on the separation distance 
R and time t, i.e., $(f, R, r) = ^{t,R). Clearly, a large 
increase of $(i, R) above the level of can lead to a 
strong grows in the frequency of the particle collisions. 

In the analytical treatment of the problem we use the 
standard equation for n{t, r): 



dn{t, r) 
dt 



W ■ [n{t,r)v{t,r)] = DAn{t,r), (1) 



where D is the coefficient of molecular (Brownian) dif- 
fusion. We study the case of small yet finite molecular 
diffusion D of particles. The equation for 6(t, r) follows 
from Eq. (1): 



deit,r) 
dt 



+ [v(i,r) • V]e(t,r) = -e(i,r)divv(t,r) 



+DAe{t,r) . 



(2) 



To study the clustering instability we use Eq. (2) with- 
out the source term oc — A'^divv, describing the effect 
of an external source of fiuctuations. Particle cluster- 
ing can also occur due to this source of fluctuations of 
particle number density. Such fluctuations were studied 
in [27, 28, 45]. In the present study we considered the 
particle clustering due to the clustering instability. Parti- 
cle clustering caused by the self-excitation of fluctuations 
of particle number density (the clustering instability) is 
much stronger than that due to the source of fluctuations 
of particle number density. 

One can use Eq. (2) to derive equation for ^{t,R) by 
averaging the equation for Q{t, r -I- R)0(i, r) over statis- 
tics of the turbulent velocity field v(i, r). In general this 
procedure is quite involved even for simple models of the 
advecting velocity fields (see, e.g., [35]). Nevertheless, 
the qualitative understanding of the underlying physics 
of the clustering instability, leading to both, the exponen- 
tial growth of <I>(t, R) and its nonlinear saturation, can 
be elucidated by a more simple and transparent analysis. 

Let us consider turbulent flow with large Reynolds 
numbers, TZe 1. Therefore, the characteristic scale 
L of energy injection (outer scale) is much larger than 
the length of the dissipation scales {viscous scale rf) 
L ^ rj. In the so-called inertial interval of scales, where 
L > r > r], the statistics of turbulence within the Kol- 
mogorov theory is governed by the only dimensional pa- 
rameter, e, the rate of the turbulent energy dissipation. 
Then, the velocity u(r) of turbulent motion at the charac- 
teristic scale r (referred below as r-eddies) may be found 
by the dimensional reasoning: u{r) « {srY^^ (see, e.g., 
[1, 46, 47]). Similarly, the turnover time of r-eddies, r(r), 
which is of the order of their life time, may be estimated 
as r(r) w r/u{r) « £~^/^r^/^. 

To elucidate the clustering instability let us consider a 
cluster of particles with a characteristic scale £ moving 
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with the velocity Vci(t). The scale ^ is a parameter which 
governs the growth rate of the clustering instability, 7. 
It sets the bounds for two distinct intervals of scales: 
L > r > £ and £ > r > r/. Note also that we cannot 
consider scales which are smaller than the size of parti- 
cles. Large r-eddies with r > £ sweep the ^-cluster as a 
whole and determine the value of Vci(t). This results in 
the diffusion of the clusters, and eventually affects their 
distribution in a turbulent flow. 

On the other hand, the particles inside the turbulent 
eddies are carried out to the boundary regions between 
the eddies by the inertial forces. This mechanism of the 
preferential concentration acts in all scales of turbulence, 
increasing toward small scales. The role of small eddies 
is multi-fold. First, they lead to the turbulent diffusion 
of the particles within the scale of a cluster size. Sec- 
ond, due to the particle inertia they tend to accumulate 
particles in the regions with small vorticity, which leads 
to the preferential concentration of the particles. Third, 
the particle inertia also causes a transport of fluctuations 
of particle number density from smaller scales to larger 
scales, i.e., in regions with larger turbulent diffusion. The 
latter can decrease the growth rate of the clustering in- 
stability. Therefore, the clustering is determined by the 
competition between these three processes. 

Let us introduce a dimensionlcss parameter (Tv, a de- 
gree of compressibility of the velocity field of particles, 
v{t, r), defined by 



_ ([divv]2) 



V X V|2) 



(3) 



This parameter may be of the order of 1 (see [32]). One 
of the reasons for the clustering instability is the parti- 
cle inertia which results in the parameter CTv ^ 0. The 
particle response time is given by 



m„ 



Gttu pa dpv 



(4) 



and the particle mass rUp is nip = (47r/3) pp. The ratio 
of the inertial time scale of the particles (the Stokes time 

scale Tp) and the turnover time of ry-cddies in the Kol- 
mogorov micro-scale T^rj) — ri/u{rf) = /v, is of primary 
importance, where u{jii) is the characteristic velocity of 
77-scale eddies. The ratio of the time-scales Tp and r(ry) 
is the Stokes number: 



St 



Tp _ 2ppa'^ 
T{rj) 9/9 77^ 



(5) 



For Tp <C r(ry) all particles are almost fully involved in 
turbulent motion, and one concludes that u{t, r) Ri v{t, r) 
and v{£) « u{£). The compressibility parameter a^, of 
particle velocity field for St ^ 1 is given by: 



2pp 
%p 



(6) 



(see [32, 35]). For small Stokes number, the clustering 
instability has been investigated in [35]. The character- 
istic scale of the most unstable clusters of small particles 
is of the order of Kolmogorov micro-scale of turbulence, 
77. The characteristic growth rate of the clustering in- 
stability is of the order of the turnover frequency of ry- 
eddies, 1/t(?7) (see [35]). In the present study we extend 
the theory of particle clustering instability to the case 
St > 1, i.e., when the particle Stokes time is larger than 
the Kolmogorov time scale, but is much smaller than the 
correlation time at the integral scale of turbulence. We 
may expect that for St > 1 the compressibility parameter 
(Tv of particle velocity field is given by: 



St' 



iSt^ 



(7) 



where a ~ 1. 



III. THE PARTICLE VELOCITY FIELD FOR 

St>l 



The equation of motion of a particle reads: 

^ = l[u(,r)-v(,r)], 
at Tp 



(8) 



where the total time derivative (d/dt) takes into account 
the time dependence of the particle coordinate r: 



di 



+ v(t, r) • V 



Now Eq. (8) takes the form: 



|+v(^,r).V 



lU(t,r) =u(i,r) 



(9) 



(10) 



In the following we analyze this equation for particles 
with the time Tp which is larger than the turnover time of 
the smallest eddies in the Kolmogorov micro-scale t(?7), 
but is smaller than the turnover time of the largest eddies 
t{L). Denote by the characteristic scale of eddies for 
which 



Tp — t{^£^ 



(11) 



This scale as well as the particle cluster scale was intro- 
duced in [30]. Note that £^/r] = St^/^. The eddies with 
£ ^ £» almost fully involve particles in their motions, 
while the eddies with £ <^ £^, do not affect the particle 
motions in the zero order approximation with respect to 
the ratio [t{£)/tp] <C 1. Therefore it is conceivable to 
suggest that the main contribution to the particle veloc- 
ity is due to the eddies with the scale of £ [which we 
denote as ve{t,r)] that is of the order of £, and much 
larger then the Kolmogorov micro-scale. Velocity ve{t, r) 
cannot be found on the basis of simple dimensional rea- 
soning because the problem at hand involves a number of 
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dimensionlcss parameters like C^/ri, etc. The main 
difficultly in determining this velocity is that in this case 
one has to take into account for a modification of the 
particle response time Tp by the turbulent fluctuations. 
The physical reason for that is quite obvious: the time Tp 
is determined by molecular viscosity of the carrier fluid 
while the main dissipative effect for motions with £ > r; is 
due to the effective "turbulent" viscosity. In order to de- 
termine the velocity V£{t,r) wc can use the perturbation 
approach to Eq. (10) (see, e.g., [48, 49]). The details of 
this derivations are given in Appendix A. This analysis 
yields the scalings of the particle velocity for St 1: 



10/9 



-{£) 



-,5/3 



(12) 



IV. THE CLUSTERING INSTABILITY OF THE 
SECOND MOMENT OF PARTICLE NUMBER 
DENSITY 

In this section we will perform a quantitative analysis 

for the clustering instability of the second moment of 
particle number density. To determine the growth rate 
of the clustering instability let us consider the equation 
for the two-point correlation function $(t, R) of particle 
number density: 

^ = [B{R) + 2U{R) ■ V + Do,0{R)VaV0] ^t, R)(,13) 

(see [35]). The meaning of the coefficients B{R), U{R) 

and Dai3{R) is as follows (for details see Appendix B). 
The function B{R) is determined by the compressibility 
of the particle velocity held and it causes the generation 
of fluctuations of the number density of particles. The 
vector U{R) determines a scale-dependent drift veloc- 
ity which describes a transport of fluctuations of particle 
number density from smaller scales to larger scales, i.e., 
in the regions with larger turbulent diffusion. The latter 
can decrease the growth rate of the clustering instability. 
Note that \5{R = 0) = whereas B{R = 0) 9^ 0. For in- 
compressible velocity field U(Jl) = and B{R) = 0. The 
scale-dependent tensor of turbulent diffusion Daj3{R) is 
also affected by the compressibility. In very small scales 
this tensor is equal to the tensor of the molecular (Brown- 
ian) diffusion, while in the vicinity of the maximum scale 
of turbulent motions this tensor coincides with the regu- 
lar tensor of turbulent diffusion. 

Thus, the clustering instability is determined by the 
competition between these three processes. The form of 
the coefficients B{R), U{R) and DapiR) depends on the 
model of turbulent velocity field. For instance, for the 
random velocity with Gaussian statistics of the particle 
trajectories these coefficients are given in Appendix B. 

Let us study the clustering instability. We consider 
particles with the size r?/\/Sc <SC a ^r], where Sc = u/D 
is the Schmidt number. For small inertial particles ad- 
vected by air flow Sc 1. There are three characteristic 



ranges of scales, where the form of the solution of Eq. (13) 
for the two-point correlation function ^(t, R) of the par- 
ticle number density is different. These ranges of scales 
are the following: (i) the dissipative range a < £ < r], 
where the molecular diffusion term oc 1/Sc is negligible; 
(ii) the first part of the inertial range rj < £ < £^ and (iii) 
the second part of the inertial range £^, ^ £ ^ L, where 
the functions B{R) and U{R) are negligibly small. 

Consider a solution of Eq. (13) in the vicinity of the 
thresholds of the excitation of the clustering instability. 
The asymptotic solution of the equation for the two-point 
correlation function <I>(t, R) of the particle number den- 
sity is obtained in Appendix B. In the range of scales 
a < £ < rj, the correlation function ^{t,R) in a non- 
dimensional form reads 



$(i?) = AiR-^^ sin(^d| In R\ + (pa) ., 



(14) 



and in the range of scales 77 < £ < it is given by 

^R) = A2R~^sm{fi\n R + (p) , (15) 

where the parameters and iid are given by Eq. (B.8) 
and the parameters A and 11 are given by Eq. (B.13) 
in Appendix B. Here R is measured in the units of rj 
and time t is measured in the units of r,, = t{£ = rj). 
We have taken into account that the correlation function 
$(i?) has a global maximum at f? = a, i.e. the normal- 
ized correlation function of the particle number density 
^{t, R = a) ~ 1. We have also taken into account that in 
the range of scales r? < £ <C £*, the relationship between 



and uj is given by: 



r{£) 



(16) 



For instance, for St ^ 1 the exponent s = 5/3 (see 
Eq. (12)). The value s = 7/4 corresponds to the 
turbulent diffusion tensor with the scaling cx R^ [see 
Eqs. (B.5)-(B.7) in Appendix B]. We consider the param- 
eter s as a phenomenological parameter. In the range of 
scales £* <C £ <C i, the correlation function $(<, R) is 
given by 

$(i?) = AsR-^' , (17) 

where A3 is given by Eq. (B.17) in Appendix B. The con- 
dition, R^^{R) dR = 0, imphes that the total number 
of particles in a closed volume is conserved. 

The growth rate of the second moment of particle num- 
ber density, the coefficients Ak and the parameters ipd^ 
P are determined by matching the correlation function 
<&(i?) and its first derivative $'(i?) at the boundary of 
the above three ranges of scales, i.e., at the points £ = rj 
and £ = £^,. For example, the growth rate 7 of the clus- 
tering instability of the second-order correlation function 
is given by 



1 



7 = 



6t^(1 -I- 3a^) 



2 (l + 3a,)2 



400 cTv 



(Ty — (Tv 
l + CTv 



(3- 



1 + 



(18) 
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where is the degree of compressibihty of the scale- 
dependent tensor of turbulent diffusion Dajji^R) (for de- 
tails, see Appendix B). Note that for the ^-correlated 
in time random Gaussian compressible velocity field, the 
coefficients B{R) and U{R) are related to the turbulent 
diffusion tensor Da0{R), i.e., 

B{R) = S/a^pDc^piR) , = Wpb^piR) , (19) 

(for details, see [32, 33, 35]). In this case the second 
moment $(t, R) can only decay, in spite of the compress- 
ibility of the velocity field. For the (5-correlated in time 
random Gaussian compressible velocity field = a^. 
For the finite correlation time of the turbulent velocity 
field 7^ a^. and the relationships (19) arc not valid. 
The clustering instability depends on the ratio a^/uy. 

The range of parameters (ctv, <j^) for which the cluster- 
ing instability of the second moment of particle number 
density may occur is shown in Fig. 1. The line = 
corresponds to the (5-corrclated in time random compress- 
ible velocity field for which the clustering instability can- 
not be excited. The various curves indicate results for 
different value of the parameter s. The curves for s = 7/4 
(dashed) and s = 5/3 (solid) practically coincide. The 
parameter s is considered as a phenomenological param- 
eter, and the change of this parameter from s = 7/4 to 
s = can describe a transition from one asymptotic be- 
haviour (in the range of scales ?7 < ^ < ^*) to the other 
{i* <C < L). The growth rate (18) of the clustering 
instability versus CTv for s = 5/3 and different values of 
is shown in Fig. 2. 

We have not discussed in the present study the growth 
of the high-order moments of particle number density (see 
[30, 35, 45, 50]). The growth of the negative moments 
of particles number density (possibly associated with for- 
mation of voids and cellular structures) was discussed in 
[45, 51, 52]. 

V. DISCUSSION 

Formation and evolution of particle clusters are of 
fundamental significance in many areas of environmen- 
tal sciences, physics of the atmosphere and meteorol- 
ogy [smog and fog formation, rain formation (see e.g., 
[16, 19, 20, 53-55]), planetary physics (see e.g., [56, 57]), 
transport and mixing in industrial turbulent flows, like 
spray drying and cyclone dust separation, dynamics of 
fuel droplets (see e.g., [58-60]). The analysis of the ex- 
perimental data showed that the spatial distributions of 
droplets in clouds are strongly inhomogeneous (see [18]). 
The small-scale inhomogeneities in particle distribution 
were observed also in laboratory turbulent flows (see 
[13, 17, 21]). 

In the present study we have shown that the particle 
spatial distribution in the turbulent flow fleld is unstable 

against formation of clusters with particle number den- 
sity that is much higher than the average particle number 




FIG. 1: The range of parameters ((Tv,cry) for which the clus- 
tering instability may occur. The various curves indicate 
results for s = 7/4 (dashed), s = 5/3 (solid), for s = 1 
(dashcd-dotted) and for s = 2/3 (dotted). The thin dashed 
line (Jv = Cj, corresponds to the 5-correlated in time random 
compressible velocity field. 




0.1 0.2 0.3 cj^ 

FIG. 2: The growth rate of the clustering instability versus Uv 
for s = 5/3. The various curves indicate results for different 

values of — + 5: 5 = 3 (solid), 5 = 1 (dashed), S = 0.5 
(dashed-dotted) and d = 0.3 (dotted). 



density. Obviously this exponential growth at the linear 
stage; of instability should be saturated by nonlinear ef- 
fects. A momentum coupling of particles and turbulent 
fluid is essential when the kinetic energy of fluid p(m^) 
is of the order of the particles kinetic energy mpnci(i'^), 
where \u\ ~ \v\, i.e., when mpUci ~ p. This condition 
implies that the number density of particles in the clus- 
ter rici ^ a~^(p/3pp). In the atmospheric turbulence the 
characteristic parameters are as follows: in the viscous 
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/ 2. 1/2, 
( M > /A 
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FIG. 3: Dependence of \/{'n?)/N versus parameter for 
s = 5/3 which corresponds to the solution of Eq. (13) for the 
two-point correlation function i?) of the particle number 
density obtained in Sec. IV. Various curves indicate results 
for different values of the parameter 5: 6 = 3 (solid), 6 = 1 
(dashed) and 5 = 0.5 (dashed-dotted), where a,^ = + 6. 



scale, ?7 ~ 1 mm, the correlation time of the turbulent 
velocity field isT^^(0.01 — 0.1)s, and for water droplets 
Pp/p ~ 10^ Thus, for a ~ 30 /xm we obtain rid ~ 10 
cm~^ (see [35]). Particle collisions can play also essen- 
tial role when during the life-time of a cluster the total 
number of collisions is of the order of number of par- 
ticles in the cluster. The collisions in clusters may be 
essential for rici ^ a~^{it,/a){p/3pp). In this case a mean 
separation of particles in the cluster is of the order of 
4 ~ a'^/'^Spp/l^pY^^. When, e.g., a SO^im we get 
4 ~ 5a « 150 Aim and ~ 3 x (10"* - lO'"^) cm^^. The 
mean number density of droplets in clouds N is about 10"^ 
cm~^. Therefore, the clustering instability of droplets in 
clouds can increase their concentrations in the clusters by 
the order of magnitude (see also [35]). Note that for large 
Stokes numbers the terminal fall velocity of particles can 
be much larger than the turbulent velocity. This implies 
that the sedimentation of heavy particles can suppress 
the clustering instability for large Stokes numbers. 

There is an additional restriction on the value of rid = 
\J (n^) which follows from the condition (n(t, r) n(t, r -|- 
R)) = ^ (e2)$(t,R) > 0, where $(t,R) is the nor- 
malized correlation function of the particle number den- 
sity $(t, R = a) = 1. Since the correlation function 
$(t, R) can be negative at some scale J?, this condi- 
tion implies that the maximum possible value of (0^) 
which can be achieved during the clustering instability 
is (9^)max = A^^/I^jmin. Thcrcforc, the number den- 
sity of particles in the cluster rici cannot be larger than 
n^x < Ar(l -|- l^lmin)"^^^- Using this criterion we plotted 
in Fig. 3 the dependencies Vv^max/^ versus param- 
eter av for different values of the parameter (5, where 
u^j, = (Tv -|- 5. We estimate |$|min using the solution of 
Eq. (13) for the two-point correlation function i?) of 
the particle number density obtained in Sec. IV. Note 
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FIG. 4: Dependence of {n?)/N versus parameter CTv which 
corresponds to the solution for the two-point correlation func- 
tion ^{t, R) of the particle number density for the case St < 1 
studied in [35]. Various curves indicate results for different 
values of the parameter 6: 5 = 3 (solid), 5 = 1 (dashed) and 
6 = 0.3 (dashed-dotted), where cr^j, = CTv -f- <5. The thick lines 
correspond to Sc = 10^ and thin lines correspond to Sc = 10"'. 



however, that this solution determines the linear stage 
of the clustering instability. In Fig. 3 we also take into 
account the conditions for the clustering instability. This 
condition implies that for a given parameter (t^ the clus- 
tering instability is excited when CTv > cr™'" (sec Fig. 1). 
For comparis on we also plotted in Fig. 4 the similar de- 
pendencies V^(^^max/^ versus parameter CTv using the 
solution for the two-point correlation function <!>(<, R) of 
the particle number density for the case St < 1 studied 
in [35]. 

In the present study we have considered the particle 
clustering due to the clustering instability. Generally, 
particle clustering can also occur due to the source of 
fluctuations of droplets number density / = B{R)N'^ in 
Eq. (13) for the second-order correlation function of par- 
ticle number density. This source term arises due to the 
term oc — A^divv in Eq. (2). Such fluctuations were stud- 
ied in [27, 28, 45]. 

Note that there is an alternative approach which deter- 
mines the particle clustering (see [61-63]). The particle 
number density fluctuations are generated by a multi- 
plicative random process: volume elements in the parti- 
cle flow are randomly compressed or expanded, and the 
ratio of the final density to the initial density after many 
multiples of the correlation time r can be modelled as a 
product of a large number of random factors. According 
to this picture, the particle number density fluctuations 
will be a record of the history of the flow, and may bear 
no relation to the instantaneous disposition of vortices 
when the particle number density is measured [61 63]. 
The particle number density is expected to have a log- 
normal probability distribution. When the random-flow 
model [61-63] with short correlation time is applied to 
fully-developed turbulence it predicts that the clustering 
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is strongest when St ^ 1, in agreement with numerical 
studies [19, 20]. 



VI. CONCLUSIONS 

In this study we considered formation of small-scale 
clusters of incrtial particles in a turbulent flow. The 
mechanism for particle clustering is associated with a 
small-scale instability of particle spatial distribution. 
The clustering instability is caused by a combined ef- 
fect of the particle inertia and a finite correlation time of 
the turbulent velocity field. The theory of particle clus- 
tering developed in our previous studies was extended 
to the case when the particle Stokes time is larger than 
the Kolmogorov time scale, but is much smaller than the 
correlation time at the integral scale of turbulence. We 
found the criterion for the clustering instability for this 
case. 
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APPENDIX A: VELOCITY OF INERTIAL 
PARTICLES FOR Tp > Tr, 

In order to determine the velocity V£{t,r) we can 
use the Wyld's perturbation diagrammatic approach to 
Eq. (10) in the Belinicher-L'vov representation (see, e.g., 
[48, 49]). This approach yields automatically a sensi- 
ble result allowing us to avoid an overestimation of the 
sweeping effect in an ordcr-by-ordcr perturbation analy- 
sis. However, keeping in mind that this approach is tech- 
nically quite involved, in this study we reformulated the 
derivation procedure and obtained the required results 
using a more simple procedure based on the equation of 
motion Eq. (10). 

To determine V£{t, r) we consider Eq. (10) in the frame 
moving with £-eddies, in which the surrounding fluid ve- 
locity u equals to the relative velocity of the i?-eddy at r, 
i.e., u{t, r) = Ui{t, r). Here one has to take into account 
that the i-eddy is swept out by all eddies with scales 
£' > I. At the same time the particles participate in mo- 
tions of ^'-eddies with I. Therefore, the relative 
velocity V(, of the ^-eddy and the particle is determined by 
i"-eddies with the intermediate scales, l^, < I' < £. This 
velocity is determined by the contribution of £*-eddies. 



and can be considered as a time and space independent 
constant during the life time of the £ eddy and inside 
it. Velocity in our approach is random and has the 
same statistics as the statistics of the turbulent velocities 
of ^*-eddies. Then Eq. (10) becomes 

(jp^ + 1^ r) = ue{t, r + uj) 

-Tp[vi{t,r)-W]ve{t,r) . (A.l) 

In Eq. (A.l) the velocity ug, is calculated at point r and 
the velocity vi is at i — it, t. For the sake of convenience 
we redefine here r ~ u^,t = r' as r and, respectively, 
r = r'-t-ii* t as r+u^ t. Note that Eq. (A.l) is a simplified 
version of Eq. (10) that we used in our derivations. 

1. First non-vanishing contribution to vi 

Since vi{t^ r) <C ui for £ <C ^* , we can find the first non- 
vanishing contribution to ve{t,r) in the limit [t{£) <^ Tp] 
by considering the linear version of Eq. (A.l): 

[tp^^ + 1^ Vi{t, r) = ue{t, r + u^t) . (A.2) 

In the Lo, k representation this equation takes the form: 

{iuiTp + 1) Vi{LJ, k) = Ui{ix> — k ■ Ut,,k) , (A. 3) 

that allows one to find the relationship between the sec- 
ond order correlation functions F^^{u>, k) and F^^{oj, k) 
of the velocity fields and uf. 

= Z;^^^^''^'' -k u.,k). (A.4) 

Functions F^^(w, k) and (w, k) are defined as usual, 
e.g., 

{2w)^S{lj + u') 6{k + k')F^^^{LO, k) = 
(vf{u,k)v^{ij\k')) . (A.5) 

The simultaneous correlation functions are related to 
their w-dependent counterparts via the integral / duj/2iT, 
e.g., 

Kfeik) = J ^F:f,{^,k). (A.6) 

The tensorial structure ofF"^{k) follows from the incom- 
pressibility condition and the assumption of isotropy: 

F:^,{k) = P'^^{k)Fu,e{k), (A.7) 

where P°''^{k) is the transversal projector: 

P"^(fc) = (5«/3 - k°'k'^/k^ . (A.8) 
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In the inertial range of scales the function F^^{u), k) may 
be written in the foUowing form: 



i?„"j(a;,fe) = P">'{k)FuAk)T{£)f[ojT{£)] . (A.9) 

Here the dimensionless function f{x) is normalized as 
follows: 



f{x) dx = 2^ . 



(A.10) 



Now we can average Eq. (A. 4) over the statistics of 
eddies. Denoting the mean value of some function g{x) 
as g{x) we have: 



/[(w - fe • ti*) r(£)] « 6[{(j -k-u^) T{e)] 

Here the dimensionless function f*{x) has one maximum 
at X = 0, and it is normalized according to Eq. (A. 10). 
The particular form of /* (a;) depends on the statistics of 
i?*-cddics and our qualitative analysis is not sensitive to 
this form. Thus, we may choose, for instance: 



Mx) = 2/[x^ + l] 



(A.12) 



In Eqs. (A. 11) we took into account that the characteris- 
tic Doppler frequency of ^-eddies (in the random velocity 
field u* of ^*-eddies) may be evaluated as: 



^^{e)^^{k-u.rc^u./e 



(A.13) 



This frequency is much larger than the characteristic fre- 
quency width of the function /[wt(^)] (equal to 1/t{£)), 
and therefore the function f{x) in Eq. (A. 11) may be 
approximated by the delta function 6{x). 

After averaging, Eq. (A. 4) may be written as 



P"^(fc)/,(o) e 



FuAk) ■ (A.14) 



Here we took into account that Tp ^ Iju^, that allows 
us to neglect the frequency dependence of f*{u^w/t) 
and to calculate this function at a; = 0. Together with 
Eq. (A.14) this yields 



FfAk)=P"^{k)F^4k)- 



(A.15) 



where we used the estimate /«(0) « 2, that follows from 
Eq. (A.12). 

The equation (A. 15) provides the relationship between 
the mean square relative velocity of f-separated particles, 
Vi, and the velocity of ^-eddies, ue: 



(A.16) 



2. EfFective nonlinear equation 

For a qualitative analysis of the role of the nonlinearity 
of the particle behavior in the ^-cluster we evaluate V in 
the nonlinear term, Eq. (A.l), as neglecting the spa- 
tial dependence and the vector structure. The resulting 
equation in w-representation reads: 



2^ J dLJi(kj2S{Lj + wi + u!2)Ve{u)i)Ve{uj2) , 
Vi{Lu) = J Vi{t) exp[—iujt]dt , 

ve{t) = [ Ve{u))exp[iu}t\duj . (A-.17) 
2tt J 

In the zeroth order (linear) approximation (A/^j 0) 

(A.18) 



which is the simplified version of Eq. (A. 3). This allows 
us to find in the linear approximation 



{vdt)) = j = j 2^^^2^ . (A-19) 

where Fu,({co) is the correlation functions of Ue{w): 

2nS{u + J)Fu,i{^) = {Ui{w)U(,{J)) , (A.20) 

similarly to Eq. (A. 5). 

In the limit Tp ^ t/u^, one can neglect in Eq. (A. 19) 
the w-dependence of F„,^(a;), which has the characteristic 
width (./u^ and conclude: 



.2 _ /„.2/ 



Ip- 



uj = {ufit)), 
in agreement with Eq. (A.16). 



(A.21) 



3. First nonlinear correction 

To evaluate the first nonlinear correction to Eq. (A.21) 
one has to substitute Vi{oj) from Eq. (A.18) into 
Eq. (A. 17) for AC: 



72 r 

14,1(0;) = I du)ldU2S{u)+C0l+U)2)^ 



ii^2 + 7p 



(A.22) 



Using Eq. (A.22) instead of Eq. (A.18) we obtain instead 
of Eq. (A.19) 
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(a;) = 



(u;2 + ^2)^2 



dujiduj2 <5(w + uJi + LU2) 



2tt 



In this derivation we assumed for simphcity the Gaussian 
statistics of the velocity field. This corresponds to a stan- 
dard closure procedure in theory of turbulence (sec. e.g., 
[1, 4]). Taking into account of deviations from the Gaus- 
sian statistics of the turbulent velocity field in the frame- 
work of the perturbation theory of turbulence does not 
yield qualitatively new results due to the general struc- 
ture of the series in the theory of perturbations after 
Dyson- Wyld line-resummation (see, e.g., [48, 49]). 
Now let us estimate 



4i 



L 



2/3 



(A.24) 



that is much larger than the result (A. 21) for vJq ob- 
tained in the linear approximation. This means that the 
simple iteration procedure wc used is inconsistent, since 
it involves expansion in large parameter [{(■*/ (-Y^^]- 



4. Renormalized perturbative expansion 

A similar situation with a perturbative expansion oc- 
curs in the theory of hydrodynamic turbulence, where 
a simple iteration of the nonlinear term with respect to 
the linear (viscous) term, yields the power scries expan- 
sion in TZe^ ^1. A way out, used in the theory of hy- 
drodynamic turbulence is the Dyson-W yld re-summation 
of one-eddy irreducible diagrams (for details see, e.g., 
[48, 49, 64]). This procedure corresponds to accounting 
for the nonlinear (so-called " turbulent" viscosity) instead 
of the molecular kinematic viscosity. A similar approach 
in our problem implies that we have to account for the 
self-consistent, nonlinear rcnormalization of the particle 
frequency 7^, =J> Tp{t) in Eq. (A.17) and to subtract the 
corresponding terms from AC- With these corrections, 
Eq. (A.17) reads: 

[iuj + Tp{t)]Ve{^^) = ^pUi{w) + . (A.25) 

Here AC is the nonlinear term A/L after substraction of 
the nonlinear contribution to the difference 



Ap^rp(^)-7p 



rp(^) 



(A.26) 



The latter relation actually follows from a more detailed 
perturbation diagrammatic approach. In our context it 

is sufficient to realize that in the limit Vp{t) ^ 7p one 
may evaluate Tp{t) by a simple dimensional reasoning: 



Tp{ti^ve/i, 



(A.27) 



which is consistent with Eq. (A.26). In addition, 
Eq. (A.26) has a natural limiting case Tp{t) 7^ when 



vtjt < 7p. Now using Eq. (A.25) instead of Eq. (A.18) 

we arrive at: 



(0). . ^vM^. 

^ ^ ' iu; + Tp{l) 



(A.28) 



Accordingly, instead of the estimates (A. 21) one has: 



2 2 ip 



Tp{t} M, 



2 Ip 



(A.29) 



The latter equation together with Eq. (A.27) allows to 
evaluate rj,(^) as follows: 



Hence the estimate (A.29) becomes 



4o ~ «l 



1/9 









5/3 











(A.30) 



(A.31) 



Repeating the evaluation of the nonlinear correction vf 2 
with the renormalized Eq. (A.25) we find that 



(A.32) 



This means that now the expansion parameter is of the 
order of 1, in accordance with the renormalized pertur- 
bation approach. This procedure yields Eq. (12). 



APPENDIX B: THE CLUSTERING 
INSTABILITY OF THE INERTIAL PARTICLES 

The clustering instability is determined by the equa- 
tion for the two-point correlation function $(f, R) of par- 
ticle number density (see Eq. (13)). The tensor Da0{R) 
may be written as 

D^ls{R) = 2D5cci3 + Dlp{R), 

DlpiR) = Dlf,{0) - DlpiR) . (B.l) 

The form of the coefficients B{R), U{R) and Dap{R) 
in Eq. (13) depends on the model of turbulent velocity 
field. For instance, for the random velocity with Gaus- 
sian statistics of the Wiener trajectories £,{t,r\T) these 
coefficients are given by 

B{R) « 2 / (&[0,$(t,ri|0)]&[T,C(<,r2|T)])dr(B.2) 
Jo 

U{R) « -2/ {v[0,^{t,rimb[T,^it,r2\T)])dT, 
Jo 

blp{R) « 2 / {v^[Q,i{t,r^mv0[T,i{t,r2\T)]) dT , 
Jo 

where b = divv (for more details, see [35]). Note that in 
this study we use Eulerian description. In particular, in 
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Eq. (B.2) the functions ij^fr, ^(t, r|T)] and &[t, ^(t, r|T)] 
describe the Eulerian velocity and its divergence calcu- 
lated at the Wiener trajectory (see [35, 65, 66]). The 
Wiener trajectory ^(t, r|s) (which usually is called the 
Wiener path) and the Wiener displacement p^^,{t, r\s) are 
defined as follows: 

pjt,r\s) = ^ v[T,${t,r\T)]dT + V2Dw{t-s), 

where w{t) is the Wiener random process which describes 
the Brownian motion (molecular diffusion). The Wiener 
random process w(t) is defined by the following proper- 
ties: (w(t))w=0, {wi{t + T)wj{t))v, ^ rSij, and {.. 
denotes the mathematical expectation over the statistics 
of the Wiener process. Since Va[T, ^{t, rjr)] describes the 
Eulerian velocity calculated at the Wiener trajectory, the 
Wiener displacement p^{t,r\s) can be considered as an 
Eulerian field. We calculate the divergence of the Eule- 
rian field of the Wiener displacements p^{t, r\s). 

Now we introduce the parameter which is defined 
by analogy with Eq. (3): 



V X xV 



, (B.3) 



where ea/s-y is the fully antisymmetric unit tensor. Equa- 
tions (3) and (B.3) imply that = in the case of 
^-correlated in time compressible velocity field. 

For a random incompressible velocity field with a fi- 
nite correlation time the tensor of turbulent diffusion 
DlfjiR) = T-i(ea(t,ri|0)C/3(t,r2|T)) (see [35]) and the 
degree of compressibility of this tensor is 



((Vx^^) • 



(B.4) 



Let us study the clustering instability. We consider 
particles of the size ri/y/Sc <^ a <^ t], where Sc = p/D 
is the Schmidt number. For small inertial particles ad- 

vectcd by air flow Sc ^ 1. A general form of the turbu- 
lent diffusion tensor in a dissipative range is given by 



2(2 + a,) ^rf_2(2a,- 



(B.5) 



3(l-h(T^) 



1) 



3(1 -hcr^) 



In the range of scales a < £ < r], Eq. (13) in a non- 
dimensional form reads: 



'dt 



= RH"{Cf + C|) + 2 R^'{Ud + Cf) + Bd^ ,(B.6) 



where R is measured in the units of rj, time t is mea- 
sured in the units of = t{£ = rj), and the molecular 
diffusion term oc 1/Sc is negligible. Consider a solution 
of Eq. (B.6) in the vicinity of the thresholds of the ex- 
citation of the clustering instability. Thus, the solution 
of (B.6) in this region is 



$(r) = AiR- 



(B.7) 



where Ai = ± i^a and 

2(Cf + C^) ' ' 2(Cf + C|) ' ^ 
{Cif = 4{Ba - 7) {Cf + Ci) - {Cf - + 2Udf , 



and 



CTv + 1 



[/d=(l/3)Brf. 



Since the correlation function $(i?) has a global maxi- 
mum at R = a, the coefficient Cf > C2 — 2Ud if is a 
real number (see below). 

Consider the range of scales r] < £ <^ i^. The rela- 
tionship between v'j and uj is determined by Eq. (16), 
where according to Eq. (12) the exponent s = 5/3. In 
this case the expression for the turbulent diffusion tensor 
in non-dimensional form reads 

dI^{R) = r(-^'~'^/^CiR''S^p + C2RM , (B.9) 



5 -I- 4s + 6(7^ 

9(1+(T^) 



(4s- l)(2cr^ - 1) 

9(1+(T^) 



To determine the functions B[R) and U{R) in the 
range of scales r/ < ^ ^ f * we use the general form of 
the two-point correlation function of the particle velocity 
field in this range of scales: 

{v^{t, r)vp{t + T,r + R)) = ^[^^^ - {ClR^do^p 



+ C2^i?„i?^)i?2(^-2)/3]/(r), 



CI 



(4- 



3(Tv) 



(l + s)(2av-l) 



3(l-hCTv) ' 3(l-ho-v) ■ 

Substitution this equation into Eq. (B.2) yields 

U{R) = Uo , B{R) = Bo R(^'-^y^ ,(B.10) 

where 



1 ' 



-Bo = p2 Uo 



and the coefficients f3i and /32 depend on the properties 
of turbulent velocity field. The dimensionless functions 
Bo and Uo in Eq. (B.IO) are measured in the units of 

For the (5-correlated in time random Gaussian com- 
pressible velocity field = (for details, see [32, 33, 
35]). In this case the second moment <^{t,R) can only 
decay, in spite of the compressibility of the velocity field. 
For the finite correlation time of the turbulent velocity 
field ^ (Tv and Eqs. (19) are not valid. The cluster- 
ing instability depends on the ratio a^/a^. In order to 
provide the correct asymptotic behaviour of Eq. (B.IO) 
in the limiting case of the 5-correlated in time random 
Gaussian compressible velocity field, we have to choose 
the coefficients pi and (32 in the form: 

/3i = 8(4s2 + 7s - 2)/27 , P2 = (4s + 2)/3 . 
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Note that when s < 1/4, the parameters /3i < and 
B{R) < 0. In this case there is no clustering instabihty 
of the second moment of particle number density. Thus, 
Eq. (13) in a non-dimensional form reads: 



'dt 



= R^^'-^y^[R^^"{Ci +C2) + 2 R^'{Uo + Ci) 



(B.ll) 



Consider a solution of Eq. (B.ll) in the vicinity of the 
thresholds of the excitation of the clustering instability, 
where {d^/dt)R^'^~'^^^/^ is very small. Thus, the solution 
of (B.ll) in this region is 



$(i?) = A2R- 



where A2 = A ± ifi, 



C1-C2 + 2Uo 
A = . ^ , , A* = 



C3 



2(Ci + C2) 



2(Ci + C2) 



(B.12) 



.(B.13) 



Since the correlation function $(i?) has a global maxi- 
mum at i? = a, the coefficient Ci > C2 — 2Uo if /it is a 
real number (see below). 

Consider the range of scales <^ £ <^ L. In this 
case the non-dimensional form of the turbulent diffusion 
tensor is given by 

Dlf^iR) = R-^'\CiR^5o,0 + C2RaR0) , (B.14) 
- _ 2(5 + 3aJ ~ _ 4(2a, - 1) 
' 9(l + a,) ' ' 9(l+a,)' 

and Eq. (13) reads: 



dt 



i?-2/3[i?2$"(C'^ + C2) +2m'Ci] . (B.15) 



Here we took into account that in this range of scales 
the functions B{R) and U{R) arc negligibly small. Con- 
sider a solution of Eq. (B.15) in the vicinity of the thresh- 
olds of the excitation of the clustering instability, when 
(9i>/9i)i?^/^ is very small. The solution of (B.15) is given 

by 



$(it!) = A3R- 



(B.16) 



where 



A3 = 



17- 



\Ci - C2 



Ci+C2 3 7(7^ 



(B.17) 



The growth rate of the second moment of particle num- 
ber density can be obtained by matching the correlation 
function $(i?) and its first derivative $'(i?) at the bound- 
ary of the above three ranges of scales, i.e., at the points 
i = r] and 1 = 1^,. Such matching is possible only when A2 
is a complex number, i.e., when C| > (i.e., /i is a real 
number). The latter determines the necessary condition 
for the clustering instability of particle spatial distribu- 
tion. It follows from Fig. 1 that in the range of parame- 
ters where /i is a real number, the parameter /i^ is also a 
real number. The asymptotic solution of the equation for 
the two-point correlation function (f)(t, R) of the particle 
number density in the range of scales a < ^ < is given 
by Eqs. (14)-(15). 
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